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I .   INTRODUCTION 

The  main  concern  when  approximating  either  an  ordinary  or  a  partial 
differential  equation  by  a  discrete  equation,  either  finite  difference  or 
finite  element,  is  the  degree  to  which  the  discrete  equation  solution 
approximates  the  solution  to  the  differential  equation.   The  closeness  of 
this  approximation  has  traditionally  been  considered  from  both  its 
quantitative  aspects,  such  as  relative  error  (Henrici  (1963)),  and  its 
qualitative  aspects,  such  as  behavior  of  transients,  propagation  of  fronts, 
etc.  (Haltiner  (1971)).   Analyses  based  primarily  on  consideration  of  the 
basically  qualitative  property  of  numerical  dispersion,  where  the 
discretization  process  causes  distortions  in  the  propagation  velocities 
for  different  frequencies,  are  especially  prevalent  (Arakawa  and  Lamb 
(1977)). 

In  this  paper,  we  apply  certain  ideas  derived  from  the  concepts  of 
transfer  functions  and  digital  filters  in  electrical  engineering  to  the 
study  of  the  qualitative  behavior  of  discretization  schemes.   The 
particular  model  we  shall  use  is  based  on  the  linearized,  one  dimensional 
shallow  water  equations,  without  mean  flow,  a  model  which  has  important 
application  to  the  study  of  the  process  of  geostrophic  adjustment.   (It 
will  be  evident  that  the  basic  approach,  however,  is  not  limited  to  this 
model.)   The  qualitative  effects  of  discretization  schemes  on  this 
equation  have  been  studied  by  Winninghoff  (1968),  Arakawa  and  Lamb  (1977), 
Schoenstadt  (1977),  and  others.   We  shall  show  that  use  of  the  transfer 
function  concept  leads  to  important  insights  into  the  differences  caused 
by  different  discretization  schemes  -  differences  that  are  not  fully 
evident  from  phase  propagation  considerations  only. 
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II.   TRANSFER  FUNCTIONS 

It  is  well  known  (Stremler  (1977))  that  linear,  space  (or  time) 
invariant  systems  can  be  fully  described  in  terms  of  their  so-called 
transfer  function.   (Figure  1)   That  is,  if  we  denote  the  input  to  the 
system  as  y.(x),  and  the  output  as  y  (x) ,  then  the  Fourier  transforms  of 
the  input  and  output,  denoted  respectively  as 


Y±(k)    =      j         y±(x)  e  1  X  dx 


(1) 


y  (k)  =   /    y  (x)  e  X  X  dx 
o       /     o 

J -co 

are  related  by 

yo(k)  =  $(k)y.(k)     .  (2) 

Equation  (2)  is  commonly  said  to  represent  the  system  in  the  transform 
domain.   The  equivalent  representation  in  the  physical  (time  or  space,  as 
appropriate)  domain  is 


/. 


yQ(x)  =   /    <|>(x  -  s)  y.(s)  ds    .  (3) 

In  equation  (3),  <J>(x)  is  commonly  referred  to  as  the  impulse  response  of 
the  system,  and  is  interpreted  as  the  response  of  the  system  to  an  input 
delta  function  at  x  =  0,  i.e.  <5(x). 

In  general,  <J>(k)  is  a  complex  valued  function.   The  physical  interpre- 
tation  of  both  the  magnitude  and  phase  of  <j)(k)  is  easily  arrived  at  by 
considering  the  response  of  the  system  to  a  single  sinusoidal  input  of 

arbitrary  frequency,  i.e. 

ik  x 
y.(x)  =  e  i     .  (4) 


It  is  very  straightforward  to  show  that 


ik.x 
y  (x)  =  <|>(k.)  e 

Thus,  if  we  represent  <f>(k)  in  polar  form  as 


?(k.)  =  |J(k.)|  e^     , 
equation  (5)  becomes: 

i(k.x  +  k.) 
yo(x)  =  |J(k.)|  e        X  * 

ik.(x  +  (k,/k.)) 

=  |<|>(k.)|  e  .  (5) 

Clearly,  from  (5),  we  can  interpret  the  magnitude  of  the  transfer  function 
as  the  factor  by  which  the  amplitude  of  a  sinusoid  of  a  given  frequency  is 
either  amplified  or  attenuated,  and  the  argument  of  the  transfer  function 
as  determining  the  amount  by  which  the  phase  of  the  output  sinusoid  is 
shifted  relative  to  the  phase  of  the  input  sinusoid. 

In  general,  neither  <J)(k)   nor  (k,/k)  is  constant,  and  hence,  the 
different  frequencies  in  a  given  input  are  amplified/attenuated  and  shifted 
in  phase  by  different  amounts.   Thus  the  output  signal  generally  has  its 
shape  (graph)  altered  from  the  input  signal.   (Figure  2)   This  alteration 
of  shape  is  commonly  called  distortion,  and,  based  on  our  discussion,  is 

composed  of  the  effects  of  two  actions  -  commonly  called  amplitude 

i  ^    i 
distortion,  which  is  due  to  the  deviation  of  |<f>(k)|  from  a  constant,  and 

phase  (or  delay)  distortion,  due  to  the  deviation  of  k,  from  a  linear 

function  of  k.   Clearly,  to  completely  understand  the  effect  of  a  linear 

system  on  an  input  signal,  one  must  know  both  of  these  effects. 
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In  electrical  engineering,  it  is  common  to  call  any  linear,  time 
(or  space)  invariant  device,  which  can  be  described  by  a  transfer  function 
as  discussed  above,  a  filter.   We  shall  use  the  terminology  hereafter  in 
this  paper. 

Lastly  we  note  that  the  energy  in  an  "output"  is  proportional  to 


/Jo    s"'    "**■        2tt 
00  J -co 


I         ■  •  > 

y    2(x)    dx   =  \-      /  |y      (k)|2   dk 


2w     J.c 


<f>(k)   y.    (k)T  dk 


/oo  r  oo 

|£(k)|2  dk}  {    /         ly^k)!2  dk} 
00  J-oo 


(6) 


=  { 


/OO  /         00 

|J(k)|2  dk}    /         y.2(x)    dx 
00  J  —CO 


and  thus  one  can  also  use  the  transfer  function  to  describe  the  relation 

between  the  energies  contained  in  the  "input"  and  "output"  to  a  linear, 

i^    1 2 
invariant  system.   (Note  the  quantity  |y(k)|  '  is  commonly  referred  to  as 

the  spectral  density  of  the  signal  y(x).) 
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III.   DISCRETE  FILTERS  AND  TRANSFER  FUNCTIONS 

A  great  deal  of  research  is  presently  directed  toward  the  study  of 
so-called  discrete,  or  digital  filters,  which  arise  when  the  continuous 
processes,  such  as  differentiation  and  integration,  in  an  invariant  linear 
filter  are  replaced  by  analogous  approximations  based  on  sampled,  or 
discrete  data.   Thus,  for  example, 

y(x  +  Ax)  =  y(x)  +  f(x)  (Ax)     ,  (7) 

is  a  digital  filter  approximating,  in  some  sense  (specifically  by  using 
an  Euler  forward  predictor),  the  continuous  integrator: 

£  -  f«     •  <8> 

The  transfer  functions  of  such  filters  can  be  easily  determined  using 
Fourier  transforms,  and  the  observation  that 


I 


,       v   -ikx  ,      ik(  Ax  )  %,  >.  ,ON 

y(x  +  Ax  )  e     dx  =  e        y(k)     .  (9) 


Thus,  the  transfer  function  of  the  process  given  by  (7)  is  (where  f(x)  is 

considered  the  "input"  and  y(x)  the  "output"),  after  some  algebra, 

tk    v  ik(Ax)/2 
$(k)  =  (Ax)e 


2i  sin(k(Ax)/2)   ' 

1 
(Note  the  transfer  function  of  the  continuous  process  is  <}>   (k)  =  —  .) 

IK. 

There  is,  however,  one  significant  difference  between  the  transfer 
functions  of  a  continuous  filter  and  its  digital  (discrete)  analog.   This 
difference  is  commonly  referred  to  with  the  terms  aliasing  or  folding,  and 
the  basic  work  on  it  is  commonly  attributed  to  Nyquist  (Clark  (1977)).   In 
its  simplest  interpretation,  Nyquist' s  work  implies  that  a  sampling  device, 
sampling  at  intervals  of  (  Ax  ) ,  is  incapable  of  resolving  waves  of 
frequency  greater  than  (ir/(Ax)),  and,  if  any  energy  is  actually  present 
in  a  sampled  continuous  signal  at  these  higher  frequencies,  it  will  be 
erroneously  resolved  (aliased)  into  a  frequency  lower  than  (tt/(Ax  )). 
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Thus,  while  a  continuous  filter  can,  in  principle,  react  to  any  input 
frequency,  and  hence  has  a  transfer  function  defined  for  all  k,  that  of 
a  digital  filter  is  considered  only  for  (complex)  frequencies  between 
+(tt/(Ax)),  the  so-called  Nyquist  cut-off  or  Nyquist  limit. 

Lastly,  we  note  that  if  a  digital  filter  has  a  significant  amplitude 
response  (i.e.  if  |<J)(k)|  is  not  close  to  zero)  near  the  Nyquist  limit, 
then  the  output  of  that  system  may  have  a  significant  portion  of  its 
energy  in  these  high  frequency  oscillations.   These  oscillations  often 
appear  to  the  observer  as  noise,  or  "noise-like",  and  hence,  in  discussing 
digital  filter  design,  many  books  and  articles  recommend  controlling  the 
transfer  function  response  near  this  limit.   (This  is  sometimes  referred 
to  as  "windowing".) 
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IV.   THE  SHALLOW  WATER  EQUATIONS  AS  A  FILTER 

The  basic  equations  we  propose  to  study,  using  the  ideas  discussed 
above,  are  the  linearized,  one-dimensional,  shallow  water  equations  in 

an  infinite  region,  with  no  mean  flow: 

3u  3h 

—  -  fv  +  g  —  =  0 

~  +  fu  =  0        >    (10) 

where  u  denotes  the  perturbation  velocity  in  the  x  direction,  v   the 
perturbation  velocity  normal  to  the  x  direction,  H  and  h   the  mean  and 
perturbed  heights  of  the  free  surface,  respectively,  and  g  >  0  and   f  >  0 
the  gravitational  and  Coriolis  parameters,  respectively.   This  model  is 
especially  important  in  the  study  of  the  meteorological  process  called  geo- 
strophic  adjustment,  which  has  been  studied  in  some  detail  from  several 
approaches  by  Rossby  (1937),  Cahn  (1945),  Winninghoff  (1968),  Blumen  (1972), 
Arakawa  and  Lamb  (1977)  and  Schoenstadt  (1977) .   This  process  of  geostrophic 
adjustment  is  important  because  it  is  the  primary  mechanism  by  which  an 
atmosphere,  modeled  by  the  so-called  meteorological  primitive  equations, 
reacts  to  errors  (either  numerical  or  due  to  errors  in  observational  data) 
in  the  initial  data.   The  dispersive  wave  nature  of  (10)  is  the  primary 
mechanism  by  which  these  errors  (commonly  referred  to  as  "imbalances")  event- 
ually spread  out  over  the  domain  of  interest,  until  the  system  reaches  a  so- 
called  "balanced"  or  "geostrophically  adjusted"  condition.   Forecasts  at  times 
before  this  balance  is  reached  are  generally  inaccurate.   Thus  a  primary  con- 
sideration in  numerically  modeling  (10)  is  to  insure  that  the  numerical 
balance  mechanism  both  fairly  accurately  models  the  physical  balancing  pro- 
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cedure,  and  that  it  proceeds  relatively  quickly. 

As  discussed  in  Schoenstadt  (1977) ,  equations  (10)  can  be  easily 
solved  by  using  spatial  Fourier  transform.  If  we  continue  to  denote 
Fourier  transforms  by  an  overhead  tilde,  e.g.: 

oo 

-ikx 


u(k,t)  =  J  u(x,t)e     dx  , 


(11) 


etc.   Then  (10)  reduces  directly  to: 

du 


—  =  f v  -  ikg  h  , 

dv    c   ~ 

dT  =  "f  u  ' 

dh     n  ,t  ~ 

—  =  -xk  H  u  , 


(12) 


together  with  initial  conditions: 


u.  =  u(k,0)  =  /  u(x,0)e 


-ikx 


dx 


(13) 


etc.   Then  since  (12)  represents  simply  a  set  of  coupled,  constant  coefficient 
(in  the  sense  of  t)  ordinary  differential  equations,  their  solution  can  be 
easily  shown,  by  the  usual  eigenvalue/eigenvector  approach,  to  be: 


f  v, 
u(k,t)  =  u.  cos  vt  H 


ikg  h. 


sin  vt  - 


sin  vt  , 


v(k,t)  ;    77  u..  sin  vt  +  { — § h  —  cos  vt}  v4  +  — §—  S1  _  cos  vt>  h. 


f 

v   1 


h(k,t)  =  - 


where : 


ikH 


u.  sin  vt  - 

l 


v 
ikHf 


v 


< 1  -  COS  Vt  >  V. 


ikgf 
2 

~2  +  — 2  C°S    I  i 


v2  =  f2  +  k2gH  =  f2(l  +  A2k2),  X  =  /ilf/f  . 
We  now  consider  (14)  in  view  of  our  discussion  of  filters  and  transfer  func- 
tions.  Observe  that  if  we  denote 


yi  = 


v. 

i 


and 


y0  = 


(14) 


(15) 


(16) 
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respectively,  as  the  "input"  and  "output"  to  the  "system"  described  by 
(10),  then  we  can  express  (14)  as 

?  =  £(k,t)y,  (17) 


where, 


<t>(k,t)  = 


cos  vt 

f 


V 


ikH 


sin  vt 


sin  vt 


v 


sin  vt 


k2gH   f2 

— S—  +  — p  cos  vt 

v     v 
ikHf 


{1  -  cos  vt} 


ikg^ 


ikgf 


sin  vt 


{1  -  cos  vt} 


v 


V 


2     2 
f    k  gH 
— 2  +  — 2   cos  vt 

v     v 


Immediate  comparison  of  (17)-(18)  with  (2)  is  not  possible  since  (J>(k,t) 
depends  not  only  on  the  transform  variable  ("k"),  but  also  on  one  of  the 
physical  variables  ("t").   However,  since  sin  vt  and  cos  vt   can  be 
expressed  in  terms  of   e      ,  and  since  k  is  real,  it  is  easily  seen 
that  the  sinusoidal  terms  in  (18)  express  only  the  phase  relationships 
(which  clearly  change  with  time),  while  the  other,  time  independent, 
coefficients  reflect  the  amplitude  effects  of  the  "filter".   In  fact,  it  is 
fairly  easy  to  show  (Schoenstadt  (1977))  that  the  sinusoidal  terms  in  (18) 
produce  dispersive,  transient  waves,  propagating  throughout  the  region. 
(These  waves  are  clearly  dispersive  since   (v(k)/k)   is  not  constant.)   The 
time  independent  coefficients  in  (18)  still  function  as  amplitude  distortion 
terms  in  the  sense  that  they  redistribute  the  energy  that  is  in  the  input 
(initial  condition)  waves,  in  a  manner  that  varies  with  wave  number,  into  the 
output  (solution)  waves. 

A  close  examination  of  (18)  shows  that,  in  fact,  the  "amplitude  distor- 
tion" in  this  "system"  is  really  governed  in  each  term  by  precisely  one  of 
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the  three  factors, 

1  ,  k  k 

v  v        °r       ~2    ' 

v 

or  the  square  of  one  of  these.   The  amplitude  response  curves  for  these  three 
expressions  are  shown  in  Figure  3.   It  is  clear  that  the  terms  with  the  co- 
efficient (1/v)   have  their  low  frequencies  least  affected,  the  terms  with 

coefficient   (k/v)   their  high  frequencies  least  affected,  and  those  with 

2 
coefficient   (k/v  )   their  frequencies  in  some  middle  range  least  effected. 

Such  responses  are  often  referred,  respectively,  to  as  low  pass,  high  pass, 

and  band  pass  filters. 

Since  the  "system"  described  by  (10)  involves  both  spatial  and  temporal 
variation,  it  is  well  known  that  not  only  the  phase  velocity   (v(k)/k)  ,  but 
also  the  group  velocity   ("Tj?)   is  important  in  describing  the  behavior,  since 
it  is  the  latter  which  governs  the  rate  at  which  energy  actually  propagates. 
Thus,  in  Figure  4,  we  display  both  of  these  velocities  for  the  system  (10). 

We  note  that  additional  insight  into  some  of  the  qualitative  behavior 
of  the  solution  to  (10)  can  be  obtained  by  decomposing  (18)  into  its  transient 
and  steady  state  components,  i.e.,  writing 


y0  -  *s  y±  +  <|>T  y.  ,       (19) 

where  <f)   and  $        can  be  obtained  by  inspection,  since  only  time  dependent 

terms  contribute  to  the  transients,  however,  we  shall  not  pursue  this  analysis 

here. 

As  a  final  comment,  note  that  the  total  energy  of  the  system  represented 

by  (10)  is  given  by: 

u2  +  v2  +  |  h2  ( 

Thus  the  spectral  density  of  the  steady  state  field  can  be  shown  to  be 

2  2 

I   |2.i   |2,gi,  i2   gHk   I   1 2   gf   i,  i2  /on>. 

IUSI   +  'Vsl   +flhSl       2   lVil   +f-2lhil   '         <20> 

v  Hv 
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Since  the  term  (k/v)  acts  as  a  high-pass  filter,  and  the  term  (f/v)  as  a 
low  pass-filter,  we  see  that  the  energy  in  the  steady  state  field  is  deter- 
mined by  a  combination  depending  primarily  on  the  high  frequency  energy  in 
the  initial  v  field,  and  the  low  frequency  energy  in  the  initial  h   field, 

We  shall  not  proceed  any  further  with  this  analysis,  since  our  main 
interest  is  to  analyze  how  discrete  schemes,  and  the  corresponding  "filters" 
they  produce,  approximate  the  continuous  (differential)  model.   However, 
there  are  certainly  other  insights  into  the  general  behavior  of  the  physical 
model  still  possible  from  this  filter  point  of  view. 
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V.   TRANSFER  FUNCTIONS  FOR  QUASI- DISCRETE  FINITE  DIFFERENCE  AND  FINITE 
ELEMENT  METHODS 


In  this  section,  we  extend  the  analysis  of  our  continuous  model  to 

the  quasi-discrete  case,  where  continuous  spatial  operations  are  replaced 

by  operations  only  at  sample  points,  while  the  continuous  time  operations 

are  retained.   This  approach  follows  that  in  Winninghoff  (1968)  and  Arakawa 

and  Lamb  (1977).   However,  we  extend  our  analysis  to  include  a  more  general 

case  of  quasi-discrete  schemes  than  they  considered.   Specifically,  we 

define  an  invariant  linear  discretization  operator  about  the  point  x  as 

any  operator  having  the  form 

n2 
L[f(x)]  =  I        c  f(x  +  n(Ax))  ,  (21) 


n 
n=-n 


where  the  c   denote  real  constants, 
n 


The  discretization  schemes  we  shall  investigate  are  specifically  those 
satisfying  the  following  conditions: 

(a)  All  discretization  operators  are  linear  and  time  and  space  invariant, 

(b)  All  discretization  operators  operate  only  on  symmetric  sets  of 
sample  points,  i.e.  in  (21),  n.  =  n_  . 

(c)  The  discretization  operators  are  balanced,  i.e.  in  (21)   c   =   c 

n      -n 

(d)  The  operators  are  applied  to  (8)-(10)  consistently,  i.e.  if  one  of 
the  derivative  terms  on  the  right  hand  side  of  (8)-(10)  is  replaced 
by  a  certain  difference,  then  all  derivative  terms  on  the  right  hand 
side  of  (8) -(10)  are  replaced  by  the  same  difference  operator. 

(e)  The  discretization  must  be  exact,  i.e.  produce  exactly  the  same 
result  as  the  analogous  continuous  operation,  on  all  at  least  arbi- 
trary linear  functions,  i.e.   f(x)  =  ax  +  b,   a,  b,  arbitrary. 
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With  these  restrictions,  it  is  fairly  simple  to  show  that  the 
quasi-discrete  schemes  we  wish  to  consider  can  all  be  reduced  to  the 
general  form: 

|^  (La[u(x,t)]}  =  f  Lg[v(x,t)]  -  g  La[h(x,t)] 

!■£  {LjvCx.t)]}  -  -f  LB[u(x,t)] 

|-  {L  [h(x,t)]}  =  -H  L  [u(x,t)] 

at   a  a 

where  L  |  ],  L.[  ],  and  L     ,  are  invariant  linear  discretization  operators 
a      0         a  r 

having  the  following  general  form: 

a 

L  [f(x)]  =  a  f(x)  +  I   a  [f(x  +  nAx)  +  f(x  -  nAx) ] 

n=l 
n„ 


L  [f(x)]  =  0of(x)  +  I   0n[f(x  +  nAx)  +  f(x  -  nAx)]    J>    (23) 

n=l 
n 

1    ° 
L  [f(x)]  =  -i  I   an[f(x  +  nAx)  -  f(x  -  nAx)] 


Ax 

n=l 

The  various  values  of  a  ,  3  ,  and  a   will  be  referred  to  as  the  filter 

n   n        n 

weights.   Thus,  for  example,  the  finite  difference  mesh  which  Arakawa  and 
Lamb  (1977)  call  Scheme  C  (Figure  5) ,  and  which  is  given  in  finite  difference 

form  as 

_3u  _  v(x+d/2,t)  +  v(x-d/2,t)  _  h(x+d/2,t)  -  h(x-d/2,t) 
dt  2  8  d 

^v  =  _  u(x+d/2,t)  4-  v(x-d/2,t) 
3t  2 

9h      u(x+d/2,t)  -  u(x-d/2,t) 
__  _H  _ 

corresponds  to  the  following  filter  weights: 

aQ  =  1,  3,  =  -^,    o     =   -r-  ,  all  others  zero.   (Note  that  in  this 
formulation  Ax  =  d/2  .) 

Determination  of  the  filter  properties  of  the  system  represented  by  (22) 
requires  taking  the  Fourier  transform  (in  x)  of  those  equations.   However,  the 
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where 


shift  property  of  these  transforms  (equation  (9))  simplifies  the  terms 
involving  the  discretization  operators,  since,  using  (9) 

{L[f(x)]}  =  \l        cn  eikn(Ax)  f (k)  .         (24) 

Thus,  the  Fourier  transform  of  (22)  is 

a(k)^-  =  f  g(k)v  -  i  g  a(k)h  , 

a(k)^  =  -f  B(k)u  , 

a(k)^r  =  -i  Ha(k)u  , 

n 

a 

a(k)  =  a  +  2  7  a   cos(nkAx)  ,  (26) 

o     *%  n 
n=l 

n3 
3(k)  =  $0  +  2  I   3n  cos(nkAx)  ,  (27) 

n=l 

n 

2    a 
a(k)  =  — -  J   a  sin(nkAx)  .  (28) 

Ax  L-    n 
n=l 

Equations  (25)  can  be  solved,  using  the  same  eigenvalue/eigenvector  approach 

as  was  used  for  (12),  to  yield: 

~   ~      a    f3%     A    i°gh0 

u  =  u^  cos  vt  H sm  vt sin  vt 

o  V  V 

f3u„  2  T7   ^2n2  .  cn 

v  = sin  vt  +  { — *—  +  — —   cos  vt}v  +  — ^-x—   {1  -  cos  vt}li        >     (29) 

v  2      2  o      2  o 

v      v  v 

~     igHuo   .   a^   iaHf_3n       A  ,~    (fV  ,a^gH     a)~ 
h  =  -  sin  vt tH.1  -  cos  vt}v  +  < — ^—  H °-^   cos  vt>ri 

v  I  v       v        I 


where  now 


a2  .  f  3  (k)  +  gHo  (k) 

V  2 

a  (k) 
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or 


v  =  ^y  /g2(k)  +  X2o2(k),  (30) 

=  v(]0 
o(k)  ' 


where 


-fA2^-2-2- 


v(k)  =  f/3  (k)  +  A  a  (k)  .  (31) 

Comparison  of  (29)-(31)  with  (14)-(15)  show  that  the  main  impacts 
of  the  quasi-discretization  described  by  (22)  are: 

(1)  The  filter  coefficients 

1   k   _k 

v  '  v  '   2 

v 

and  the  squares  of  the  first  two,  are  replaced,  respectively, 

8  a       go 

v  '  v  '   2 
v 

and  the  squares  of  the  first  two  of  these  expressions. 

(2)  The  phase  velocity 

v   f 


is  replaced  by 


■    A  +  A2k2 


f       /2~TT 

i — 7TT  ^3      +  A   q 
ka(k) 


Since  the  filter  coefficients  discussed  above  determine  how  the 
energy  in  the  initial  disturbance  is  partitioned  into  the  transient  and  steady 
state  fields,  then  the  degree  to  which  the  discrete  filter  coefficients  approxi- 
mate the  continuous  ones  can  be  viewed  as  a  measure  of  the  degree  to  which  the 
discretization  accurately  protrays  the  continuous  model  energy  distribution. 
Any  difference  between  a  continuous  (differential)  coefficient,  and  the  corres- 
ponding discrete  coefficient,  should  be  considered  as  introducing  a  distor- 
tion between  the  continuous  and  discrete  solutions. 
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Similarly,  any  differences  between  the  continuous  phase  and  group  velocities 
and  the  discrete  ones  will  result  in  distortion  due  to  spreading  of  the  waves, 
This  effect  will  be  especially  noticeable  if  the  distortion  is  pronounced  at 
(spatial)  frequencies  that  carry  significant  energy,  i.e.  for  which  filter 
coefficients  are  "large". 

Actually,  as  noted  before,  the  phase  velocity  is  less  important  than 
the  group  velocity, 


in  the  continuous  case,  and 


dv 
dk 


dv 


dk 
in  the  discrete  case,  since  this  quantity  measures  the  velocity  at  which  the 
energy  in  the  different  frequencies  propagates. 
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VI.   COMPARISON  OF  THE  ARAKAWA  SCHEMES 

The  specific  discretization  schemes  which  we  shall  analyze  using 
(22)-(30)  are  based  on  the  arrangements  of  sample  points  described  in 
Arakawa  and  Lamb  (1977)  as  Schemes  A-D.   The  arrangements  of  these  points 
in  one  dimension  is  shown  is  Figure  5.   We  shall  consider  both  finite 
difference  and  finite  element  formulations  of  the  discrete  problem, 
whereas  Arakawa  and  Lamb  considered  only  second  order  finite  differences. 
The  finite  difference  formulations  are  derived  for  both  second  and  fourth 
order  approximations  to  the  spatial  derivative  terms,  while  the  finite 
element  formulations  are  derived  (for  Schemes  A-C  only)  using  piecewise 
linear  elements  (the  so-called  "hat"  function  basis).   The  filter  weights 
derived  for  each  of  these  methods  are  shown  in  Tables  1-3.   Note  that  since 
each  scheme  "samples"  the  initial  disturbances  at  a  distance  of  d,  then  the 
shortest  wave  that  can  be  resolved  in  each  case  is  of  frequency  (u/d) . 

In  Figures  6-8,  the  amplitude  coefficients  of  the  discretized  approxi- 
mate "filters"  are  compared  to  the  "filters"  in  the  differential  model,  for 
the  values  of  the  physical  parameters  given  in  Arakawa  and  Lamb  (1977) .  The 
following  qualitative  generalizations  are  clear: 

(1)  In  each  of  the  methods  (second  order  difference,  fourth  order 
difference,  or  finite  elements),  the  arrangement  of  points  specified  by 
Scheme  A  appears  to  produce  the  greatest  distortion,  and,  furthermore,  this 
distortion  is  most  pronounced  at  high  frequencies  (short  waves). 

(2)  In  each  of  the  methods,  Schemes  B,  C  and  D  very  accurately  approx- 
imate the  transfer  function  for  the  "high  pass"  filter. 

(3)  In  each  of  the  methods,  Scheme  C  understates  the  amount  of  energy 
distributed  into  the  short  waves  by  both  the  low-pass  and  band  pass  filters, 
while  Scheme  B  overstates  this. 
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(4)  In  all  cases,  the  finite  element  formulation  appears  more 
"accurate"  than  fourth  order  differences,  which  in  turn  is  more  "accurate" 
than  the  second  order  differences. 

(5)  The  (group)  velocity  distortion  is  especially  serious  in  Schemes 

A  and  D,  since  for  short  waves  there  is  actually  a  reversal  in  the  direction 
of  propagation. 

Note  that  even  though  Scheme  B  seems  to  have  the  least  "distortion"  in 
each  case,  there  is  still  significant  phase  and  group  velocity  distortion 
near  the  cut-off.   Thus,  we  should  not  conclude  immediately  that  Scheme  B 
will  produce  the  least  "noisy"  solutions.   This  is  because  the  apparent  noise 
in  the  solution  is  a  combination  of  both  amplitude  and  phase  behavior,  and 
any  phase  distortion  may  lead  to  perceived  "noise"  in  the  solution  if  (as 
is  the  case  in  Scheme  B)  significant  energy  is  retained  (due  to  the  amplitude 
distortion  effects)  in  the  short  waves.   In  fact,  it  may  be  preferable  to 
tolerate  more  amplitude  distortion,  if  it  acts  to  reduce  the  amount  of 
energy  in  the  short  waves,  which  are  generally  the  slowest  propagating,  and 
which  produce  the  appearance  of  noise  in  the  solutions. 

As  a  final  check  on  our  analysis,  we  compared  the  height  fields  for 
both  the  continuous  and  all  the  discrete  models,  using  the  same  methodology 
as  in  Arakawa  and  Lamb  (1977).   An  initial  disturbance  given  by: 

ui(x,0)  =  0 


h  (x,0)  =  0 


v.(x,0)  =  . 


vn   ,  -  a  <  x  <  a 


0    ,  otherwise 

was  considered,  with  the  physical  parameters  having  the  following  values 

-2 
g  =  10  m  sec 

H  =  103  m 

f  =  10~4  sec"1 

a  =  d  =  A/2. 
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The  transform  of  the  initial  v  field 


2V„  sin(ka) 
~      0 
v. 


on 


i        k 
was  then  band  limited  to  the  region  -  -r  <  k  <  — ,  and  the  exact  soluti 

for  h  was  computed  as  the  inverse  transform  of 

ikHf  ._  2VQ  sin(ka) 

x—  {1  -  cos  vt}{ } 

2  k 

V 

in  the  continuous  (differential)  case,  and  the  inverse  transform  of 

ia3fH  n      ,V     2V0  sin(ka) 

r  {1  -  cos(a  t)}{ k } 

v 
in  each  of  the  discrete  cases.   In  each  instance,  the  inversion  was  accom- 

plished  using  Simpson  s  rule  on  the  interval  (0,  -r) ,  with  600  subintervals. 

The  different  h  fields  are  shown  at  Figure  9.   Note  that  these  agree 

with  our  predictions,  based  on  the  filter  analysis,  in  that: 

(1)  In  the  second  order  finite  difference  methods,  Scheme  A  displays 

a  great  deal  of  numerical  noise,  as  predicted  due  to  the  great  amplitude  dis- 
tortion of  its  filter  near  (ir/d)  .   Schemes  B  and  D  show  somewhat  less  noise 

than  Scheme  A,  but  still  more  than  Scheme  C,  this  due  to  the  manner  in  which 

2 
Scheme  C  controls  the  amplitude  near  (ir/d)  in  both  the  (a/V)  and  ($a/V  )  terms. 

(2)  Similar  comments  hold  for  the  fourth  order  finite  difference  methods, 
though,  in  general,  due  to  better  phase  propagation,  the  spreading  of  the 
energy  in  the  short  waves  is  not  so  pronounced,  and  hence  the  noise  in 
Schemes  B  and  D  is  reduced. 

(3)  In  the  finite  element  methods,  Scheme  A  still  remains  quite  noisy, 
due  to  its  amplitude  distortion  near  (ir/d)  .   Schemes  B  and  C  are  virtually 
identical,  since  Scheme  B  seems  to  compensate  for  slightly  more  energy  in 
the  short  waves  by  maintaining  better  phase  relationship. 
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VII   CONCLUSIONS 

In  this  paper  we  have  adapted  concepts  that  are  commonly  associated 
with  the  design  of  filters  in  electrical  engineering  to  the  analysis  of 
discretized  representations  of  partial  differential  equations,  and  spec- 
ifically to  the  one  dimensional  shallow  water  equations.   We  have  shown 
that  these  concepts  can  provide  valuable  insights  into  why  certain  solution 
schemes  provide  "better"  solutions  than  other,  and  a  useful  tool  for  com- 
paring both  qualitative  and  quantative  aspects  of  different  schemes. 

We  have  clearly  demonstrated  that,  in  contrast  to  the  analysis  used  in 
earlier  studies,  a  complete  analysis  of  the  suitability  of  a  certain  discrete 
scheme  to  model  a  physical  process,  in  our  case  the  process  of  geostrophic 
adjustment,  depends  not  only  on  providing  reasonably  accurate  duplication  of 
the  phase  propagation  (delay)  characteristics  of  the  modeled  system,  but 
also  on  producing  reasonably  accurate  duplication  of  the  amplitude  response 
characteristics.   Furthermore,  because  of  the  phase  distortion  near  the 
Nyquist  cut-off  of  discrete  schemes,  it  is  often  desirable  to  control  the 
amplitudes  of  short  waves. 

In  the  process  of  geostrophic  adjustment  studies,  we  have  concluded  that 
schemes  that  use  unstaggered  grid  points  are  generally  poor,  irrespective  of 
the  particular  finite  difference/finite  element  method  used,  due  to  some  in- 
herent properties  of  the  model.   The  difficulties  which  have  been  reported 
with  these  methods  in  the  literature,  and  especially  the  problem  of  "noise" 
in  finite  element  models,  appears  directly  attributable  to  the  amplitude 
response  near  the  Nyquist  cut-off  methods  based  on  an  unstaggered  arrange- 
ment of  points.   This  is  certainly  related  to  the  tendency  of  schemes  based 
on  unstaggered  grid  points  to  produce  solution  separation,  and  appears  to 
arise  from  the  fact  that  the  coupling  of  the  fields  in  the  basic  differential 
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equation  simply  makes  specifying  all  the  quantities  at  every  grid  point  an 
overspecif ication  of  the  problem. 

Lastly,  we  conclude,  that,  based  on  the  one  dimensional  evidence  only, 
the  method  identified  by  Arakawa  and  Lamb  (1977)  as  Scheme  C  has  been  quite 
popular,  not  only  due  to  its  good  phase  propagation  characteristics,  but 
also  due  to  its  tendency,  in  either  finite  difference  or  finite  element  for- 
mulations, to  "window"  out  much  of  the  high  frequency  noise  in  the  discre- 
tized  model. 

In  summary,  the  concepts  we  have  applied  yield  extremely  valuable 
insights  into  the  qualitative  behavior  of  discrete  methods  for  approximating 
partial  differential  equations.   In  many  ways,  we  have  only  started  to  bring 
much  of  the  emerging  understanding  of  digital  filters  to  bear  on  this  problem. 
A  significant,  and  interesting  outgrowth  will  be  the  consideration  of  two 
dimensional  problems. 
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TABLE  1 

FILTER  WEIGHTS  FOR  SECOND  ORDER  FINITE  DIFFERENCES 

(Ax  =  d/2) 


SCHEME      DENOTED      a       0       B,      (T,       a„ 

o       o       1       1       2 


A2 


B2 


C2 


D2 


0 

0 

1 
4 

0 

1 
2 

0 

1 

2 

1 
2 

0 

1 
2 

0 

1 

4 
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TABLE  2 

FILTER  WEIGHTS  FOR  FOURTH  ORDER  FINITE  DIFFERENCES 

(Ax  =  d/2) 


SCHEME      DENOTED      a       3       0,      a,      a0      a„     a. 

O  O  1  1  l  J         q 


A4        1       1       0       0       ^     0    --^ 


B4 


C4 


D  D4 


0 

27 
48 

0 

1 
48 

0 

1 

2 

27 
48 

0 

1 
48 

0 

1 
2 

0 

8 
24 

0 

1 

24 
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velocities  for  the  second-order  finite  difference  filters 
as  compared  to  the  continuous  (DIFF)  filter. 

Figure  7  -    (a)-(c)  Amplitude  coefficients,  (d)  phase  and  (e)  group 

velocities  for  the  fourth-order  finite  difference  filters 
as  compared  to  the  continuous  (DIFF)  filter. 
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